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1. Introduction 

A density operator of a quantum system is determined by any informationally complete set 
of measurements performed on the system; this means that the state is uniquely specified by 
the collective outcome statistics of such measurements (see, for instance [22], [6]). However, this 
point of view is rather abstract; in practical applications one instead aims to derive explicit 
reconstruction formulas for the density operator in terms of the empirical distributions in ques- 
tion. Of course, informational completeness of the measurements is necessary for the existence 
of such reconstruction formulas. In quantum optics one typically uses the set of rotated quadra- 
tures pm CGE CEB] , which can easily be measured by homodyne detection [16]. Another option is 
to use phase space observables. In particular, one can reconstruct the density matrix from the 
collection of phase space observables generated by number states. The associated distributions 
are sometimes called displaced photon number distributions, and the entire collection is called 
photon number tomogram by Manko et al. [T9l 1201 |2~T]. 

The purpose of this paper is to use the method of infinite matrix inversion, as in [H], to 
derive state reconstruction formulas involving the measurement outcome distributions of phase 
space observables generated by the number states. We consider two different types of formulas, 
involving (1) the entire tomogram, and (2) only a single observable. We find one formula of 
type (2) for each displaced photon number distribution; up to our knowledge, such formulas 
have previously been obtained only for the observable generated by the vacuum state. 

The paper is organized as follows. In section [2] we fix the notations and consider the informa- 
tional completeness of phase space observables. The relevant results concerning the inversion 
of infinite matrices are proved in section [31 Section [4] contains the main results of this paper. 
After presenting the basic properties of the phase space observables generated by the number 
states and discussing the possibility of measuring the observables, we prove two reconstruction 
formulas. In section [Sj we discuss the advantages and disadvantages of the two reconstruction 
scenarios considered in this paper. 



2. Preliminaries 



Let TC be a complex separable Hilbert space and {\n) \ n G N} be an orthonormal basis of TL 
where N := {0, 1,2,.. .}. The basis is identified with the photon number basis, or Fock basis, in 
the case where TC is associated with the single mode electromagnetic field. Let a and a* denote 
the usual raising and lowering operators associated with the above basis of TC, and let iV = a* a 
be the selfadjoint number operator. Now the phase shifting unitary operators are R(9) := e l6N . 
Define the shift operator of the complex plane D{z) = e za *~ za , z G C, for which the identities 
D(z)* = D(z)- 1 = D(-z) and R(6)D(z)R(6)* = D(ze w ) hold. The matrix elements of D(z) 
with respect to the number basis are 



(1) (m\D(z)\n) = (_i)max{o,n- m} ^ (m -n) / mm{m,n}! 2/2 |m _ re| | | 2 

' 1 1 l / max{m,n}! mm{ m ,n}\ h 



where z = re ie and 





, U\ \S-Uj 
u=0 x 



is the associated Laguerre polynomial. 

Let C(TC) be the set of bounded operators on TC, and T(TC) the set of trace class operators. 
We let || • ||i denote the trace norm of T(TC) and the operator norm of £(TC) is denoted by || ■ ||. 
When TC is associated with a quantum system, such as the single mode electromagnetic field, the 
states of the system are represented by positive operators p G T(TC) with the unit trace, density 
operators, and each state is fully characterized by the matrix elements p mn '■= ( m \p\ n ) with 
respect to the given basis. The observables are associated with the normalized positive operator 
measures (POMs) which, in the case of phase space observables, are defined on the Borel a- 
algebra £>(C) of subsets of C = M 2 Q The measurement outcome statistics of a phase space 
observable E : £>(M 2 ) — > C(Ti) in a state p are given by the probability measure X i— > tr[pE(X)]. 

For each positive operator K of trace one, define the phase space POM E K : £>(C) — > C(H) 

by 

(2) E*(X):= f D(z)KD(zy ir Z 



JX 7T 

where the integral exists in the a-weak sense. This measure is covariant in the sense that 

D{a)E K {X)D{a)* = E K {X + a), 

for all X G £>(C) and a G C. Furthermore, each covariant phase space observable is of the 
above form |T2|, [24]. We use the notation G K for the operator density related to E K , that is, 
G K : C — > £(7i), z i— > G K {z) = D(z)KD(z)* . For a fixed density operator p, the phase space 
probability measure associated with E K , i.e. X \— > tr[pE K (X)], has a density : C — > [0, oo), 
given by 

G*(z) = tr[ P D(z)KD(zT) = tr[pG K (z)}. 
If K is a one-dimensional projection, that is K = \ip)(ip\ for some ip G H, \\ifj\\ = 1, we use the 
notations E^ := E^M, G^ := G^M and Gf := df m respectively. 



1 A normalized positive operator measure, denned on a c-algebra E of subsets of a set fi, is a map E : E — > 
C(7i) which is cr-additive in the weak operator topology, and has the property E(O) = / (the identity operator), 
that is, for which X i— > tr[pE(X)] is a probability measure for each state p. 

2 



When reconstructing the state of the system directly from some measurement data, the 
measured observables are required to distinguish between any two states: 

Definition 1. A set M. of observables E : B(Q) — > C(H) is informationally complete, if any 
two states p and p' are equal whenever tr[pE(X)] = tr[p'E(X)] for all E £ M and X £ £>(f2). 

In other words, the informational completeness of a set M. of observables means that the 
totality of the corresponding measurement outcome distributions determines the state p of the 
system. Clearly, a set M. of observables is informationally complete if and only if p = whenever 
p is a selfadjoint trace class operator with tr[pE(X)] = for all E £ M. and X £ B(Q). If M. 
consists of a single observable E, we say that E is an informationally complete observable. A 
covariant phase space observable E K is known to be informationally complete if tv[KD(z)] ^ 
for almost all z £ C |2J. As a consequence of this, we get the following lemma, which shows that 
E K is informationally complete whenever if is a finite matrix. In particular, the observables 
generated by the number states are informationally complete. 

Lemma 1. Let K be a positive operator with unit trace, whose matrix representation with 
respect to the number basis {\n) \ n £ N} is finite. Then the covariant phase space observable 
E K generated by K is informationally complete. 

Proof. Since the matrix representation of K is finite, K can be written as a finite sum K = 
Ylmn=o^ mn \ m )( n \- Due to the linearity of the trace, we then have 

k 

tx[KD{z)] = K mn (n\D(z)\m) 

m,n=0 

for all z £ C. According to equation (pQ), we have (n\D(z)\m) = exactly when 

I \\m-n\r\m-n\ /i |2\ _ r, 
\ Z \ Ll mm{m,n}\\ Z \ ) ~ U ' 

which is a polynomial of \z\ of order m + n. We thus find that (n\D(z)\m) = for only a finite 
number of points z £ C, which then implies that 

tv[KD(z)\ ^ 

for almost all z £ C. Hence, E K is informationally complete. □ 

Consider now an arbitrary covariant phase space observable E K . Let P n be the projection 
onto the n-dimensional subspaces spanned by the vectors \k), k = 0, l,...,n — 1, that is 
P n = J2k=o\k)(k\- Since tv[K] = 1, there exists a smallest natural number n such that 
tr[P no KP m ] 7^ 0. For each n > n , define the truncated operator K n = tr [ P X KP \ P n KP n , where 
the normalization assures that it is a positive operator of unit trace. According to lemma CD, 
each observable E Kn is informationally complete. It is a well known fact that the sequence 
(K n ) n£ ^ converges to K in the trace norm. For each state p and X £ £>(C) we then have 

tr[pE^(X)] = / tr[pD(z)K n D(z)*}— = [ tr[K n D(z)*pD(z)} — 
Jx k Jx x 

= [ tr[K n D(z)pD( Z y]—=tr[K n EP(-X)} 

J-X 7T 

and similarily for K. This then implies that 
\tx[ P E K -{X)]-tx[pE K {X)]\ = \ti\K n E^-X)]-ti[KE^-X)}\ < \\K n - ^||i||E"(-JC)|| - 0, 
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as n — > oo. In this way, the measurement of E is obtained as a limit of measurements of 
informationally complete observables. In particular, the measurement of an informationally 
incomplete observable can be obtained as such a limit. 



3. Matrix inversion results 

In this section we prove the relevant results concerning the inversion of infinite matrices. The 
first result shows that any infinite upper triangular matrix with nonzero diagonal elements has 
a formal inverse. 

First of all, notice that the product of two or more upper triangular matrices is always a 
well defined upper triangular matrix, in the sense that the matrix elements of the product 
matrix are well defined finite sums. To clarify this, consider the matrices A = (a mn )m,neN an d 
B = (6 mn )m,neN, for which a mn = = b mn for n < m. Now the matrix elements of the product 
matrix are 

oo m+l 

(AB) 

m,m+l / &m,kbk,m+l / <^m.,fc"fc,m+Z 
fc=0 k=m 

for all I G N, and (AB) mjJl = for n < m. Similarly, any finite product of upper triangular 
matrices is well defined. If C is a strictly upper triangular matrix, that is, the diagonal elements 
are zeros, then for each m, I G N we have (C k ) m . m+ i = when k > I. In this way, the infinite 
series 

oo 
k=0 

is well defined in the sense that 

' oo \ oo I 

v fc=0 / m ,m+l k=Q k=0 

that is, the series reduces to a finite sum for each m, I G N. 

If A is an upper triangular matrix with unit diagonal elements, then the matrix (7 — A) is 
stricly upper triangular. Thus, the series 



k=0 

is well defined in the above sense. The following lemma shows, that such a series is actually 
the formal inverse of A. 

Lemma 2. Let A = ((wOm^eM be an upper triangular infinite matrix with unit diagonal, that 
is, a mn = for n < m, and a mm = 1 for all m G N, and let B = (6 mn )m,neN be a matrix for 
which b mn = for n < m and b mtm+ i = ^2 k=0 [(I — A) k ] mtin+ i for all I G N. Then A and B are 
formal inverses of each other, that is (AB) mn = 5 mn = (BA) mn . 

Proof. First notice that for n < m we have trivially (AB) mn = = (BA) mn since they involve 
empty sums. The case of the diagonal elements is also clear since for example (AB) mm = 
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a mm b mm = 1. Suppose now that n — m + I, where I > 0. Define the matrices A := (ajj)ijio 
and B := {bij)™^ as finite cut-offs of the corresponding infinite matrices. Now 

m+l 

(AB) 

m,m+l / Q>mkbk,m+l 

(AB) 

k=m 

so it is sufficient to prove the claim for finite matrices. Clearly 

m+l 



j 



k=0 



since (J - Af = when fc > m + I + 1. Thus 

m+i m+Z— 1 

AB = (J-(/-A))B = ]T(/-A) fc - 



fc+i 

fe=0 fc=0 



fe=0 fc=0 

and hence 

(AB) m)m +i = (AB) mtm +i = 

for all m G N and I > 0. In a similar fashion one proves that (BA) m ^ m+ i = for all m G N, 
/ > 0. □ 

As an immediate consequence, we find the inverse of an upper triangular matrix with nonzero 
diagonal elements. 

Corollary 1. Let A = (a mri ) mi „ G N be an upper triangular matrix with nonzero diagonal ele- 
ments, and U = (w m n)m,neN a diagonal matrix with u mm = a~^. Then A has a formal inverse 
B = (6 mn ) m ,neN such that b mn = when n < m and 

1 1 1 1 
b m , m+ i = - UA) k ] m , m+l = - AU) k U m+l 

a m+l,m+l k=Q 0"m,m fc=Q 

for all m,l G N. 

Proof. First notice that the problem again reduces to the case of finite matrices. Now UA and 
AU are upper triangular matrices with unit diagonals, so taking suitable cut-offs of these, the 
claim follows from elementary calculations. □ 

Consider now a finite sequence (c„)^ =0 C C in the above case. Define the sequence (d m ) me ^ 
via 

oo k 
n=0 n=m 

Since d m = for m > k, the sequence is actually finite. Define (4JneN via 

oo k 

c n ^ ^ b nm d m ^ ^ b nm d m . 

m=0 m=n 
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Again, c' n = for n > k, and inserting d m into the above equation gives us for n < k 

< = ± ± - ± t = t(t «~w) c = £ « = <, 



since n' < fc. This then implies that when restricted to the vector space of finite sequences, 
the linear mappings corresponding to the matrix and its formal inverse are inverse mappings 
of each other. 

The second lemma deals with a special case of an upper triangular matrix, namely, one that 
is also an infinite dimensional Toeplitz matrix. That is, for all / G N, the Ith diagonal elements 
Om,m+2, tti G N, do not depend on m. It turns out that the formal inverse (fe mri ) m ,n6N is also an 
upper triangular Toeplitz matrix. In this case we also find a sufficient condition for inverting 
the relation 

oo 
n=0 

where (c„) ne N is an infinite sequence, as 

oo 

c = b d 

m=0 

Lemma 3. Let I G N, I > 1, ao, ai, . . . , <fy G C, do 0, and define the matrix A = (a sn ) Sjn6 N 
for which a sn = a n _ s , when s < n < s + I, and a sn = otherwise. Let B = (6 ns )„ jSe N be the 
formal inverse of A. 

(a) There exist a unique sequence (b u ) ue -M C C such that b ns = b s ^ n when s > n. 

(b) Let ( c n )ngN C C and define the sequence (d s ) s6 ^ via d s — ^2^L a, sn c n . Suppose that for 
a given n G N, the condition linim-^ ak- n bm-kC m = is satisfied for k = n+1, . . . , n+l. 
Then 

oo 

Cra ^ , bnsds- 
s=0 

Proof. To prove (a), we are going to show that the matrix elements b n>n+ k, n, k G N, do not 
depend on n. According to Corollary CD, we have b ntTl+ k = — — U A) u ] nyn+ k for k > 0, 

and b ns = otherwise. First note that b nn = for all n G N. Suppose now that k > 1. Since in 
this case we have simply UA = j-A, we get (I — C/A) ns = — for s > n, and (I — LM.) ns = 
otherwise. A direct calculation now gives us 

oo oo oo 

[(L-UA)X,n+k = E E E (I-UA) ntl (I-UA) tlt2 ---(I-UA) tu _ un+k 

tl=n+l *2=ti + l t u _i=t tt _2+l 

^ \ m n+M ti+M tu-2+M 

/, /, "• /, at!-nat 2 -ti ■ ■ ■ a n+ k-t u -! 

fln / — ' * — 

7 ti=n+lt 2 =ti+l t u _i=t u _ 2 +l 
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for u > 1. After suitable changes in the summation indeces, we obtain 

_^ k , 1 \ " nJrU t u -2+U 

bn,n+k = h 2, \ ' ' ' a ti-nO>t2-ti ' ' ' a n+k-t u -! 

Cln \ Cln I 

u u =0 v u/ ti=n+lt 2 =ti+l t u _i=t u _ 2 +l 

fe /- j \ u U tl+U t u -2+U 

= — + J2 J2 '■■ a tl a t2 - tl ---a fc _ tu _ 1 , 

U =0 \ 0/ t 1 = lt 2 =ti+l t u _i=t u _ 2 +l 

which goes to show that does not depend on n. Consequently, the sequence (6/)j £ n, 

= b i is of the desired form. In addition, it is clearly unique. 

To prove (b), we first deal with the case n — 0. Consider the partial sum Sk :— J2^ =0 bn s d Sj 
for k > I — 1 (Recall the assumption I > 1.). We put in the expression 

oo s+i fe+i 

(i s = ^ ^ Cl sn 'C n r ^ ^ dsn'Gn' ^ ^ a sn'Cn'i k ^ S, 

n'=0 n'=s n'=0 

to get 

fc k+l k+l / k \ 

<Sfc = ^2 b ° s ^2 asn ' Cn ' = ^2 \ ^2 b ° sClsn ' J C n >. 
s=0 n '=0 n'=0 \s=0 / 

According to (a), the sum in parenthesis equals 5n n ', provided that the summation covers the 
full range of nonzero a sn r.s. This happens exactly when k > n' . Thus we get 

k+l / k \ k+l 



Sk = C + ^ [^2 b 0s a sn' ) C n > = C + ^ I ^2 bsCLn ' 
n'=k+l \s=0 / n'=k+l \s=n'-l 



I / k \ I / I \ 

— c o + ^ I ^ b s a n > +k _ s J c n / +fe = c + ^ I ^ 6 n / + fe_ s a s J c„/ +fe 

n'=l \s=n'+A;-i / ra'=l \s=n' / 

where the third equality is obtained by substituting n' \— > n' + k in the outer sum, and the 
fourth equality by substituting s n' + k — s in the inner sum. Suppose now that the limit 
condition holds for n = 0. Then 

= hm ^2 ( b k-s a s I c fe = Hm I fr„'+ fc _ s a s ) c n , +fc 

n'=l \s=n' / n'=l \s=n' / 

proving that lim^oo Sk = c . 

Now fix an n G N and define a translated sequence c n / — c n '_|_ n , n' G M, with 

oo s+l s+n+l 

d s := ^2 a sn'C n ' — ^2 a n'~s c n+n' = ^2 a n'-(s+n) c ri = d s+n , S G N. 
n'=0 n'=s n'=s+n 

Hence, we have the convergence 

oo oo oo oo 

Cn ^ ^ b ns d s ^ ^ b ns d s ^ ^ b n s + n d s + n ^ ^ ^os^s+n 

s=0 s=n s=0 s=0 

exactly when 

oo 

Co = ^2 b ° s ^ s 
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which, according to the result just obtained, happens if and only if lim m ^oo akb m -kC m = 0, 
k = 1, . . . , I. But this is equivalent to the claimed limit condition, and the proof is complete. □ 

Remark 1. (a) According to the proof, a necessary and sufficient condition for the conver- 
gence of the series c n = ^2^L b ns d s is that the remainder 

i i 

ri -EE 

n'=l s=n' 

goes to zero in the limit k — ► oo. This is not equivalent to the limit condition of lemma 
[3] in general. 

(b) The limit condition of lemma [3] cannot be relaxed. Indeed, at least in the case / = 1, it 
is also necessary for the convergence of the series, if we assume a\ ^ 0. This is apparent, 
since the remainder term R% contains only one term then. Another example is given 
by I = 2, ai = 0, and a 2 ^ 0. In this case, b 2s = (—^Y^i an d b 2s +i = 0. Hence 
the remainder term is K% = a 2 bk-\c n+ k+i for odd k, and = a 2 bkC n+ k+2 for even 
k. Thus, the necessary and sufficient condition for the convergence for the c n series is 
lim^oo a 2 b 2 kC2k+n+2 = 0. This is just the same as the limit condition for all sequences 
(c„)„ e N, since ai = implies that Hindoo aib 2k c 2 k +n +i = trivially. 



4. Phase space observables generated by the number states 

For each sGN, let G^ : [0, oo) x [0, 27r) — ► C(H) be the operator density associated with the 
phase space observable generated by the number state i.e. G^(r, 6) = D(re l9 )\s)(s\D(re t9 )* . 

For any state p, let G p s \r,9) := tr[pG'^(r, 9)} be the corresponding probability densit}H . The 
informational completeness of the corresponding observable 

B(C) 3 Zn E^(Z) = [ D(z)\s){s\D(z)*— e C{H) 

JZ 7T 

follows directly from lemma [Q and thus the reconstruction of the state is, in principle, possible 
from the measured distribution. 

The matrix elements of the operator density with respect to the number basis are 

(n\a s \ ri e)\m) = e^ n -^fUr), 

where 

/* m (r):= (n\D(r)\s)(s\D(ry\m). 
Thus, the probability density Gp\r,6) = tr [pG^(r, 6)} can be written as 

oo oo 

G p s \r,6)= Pmn{n\G\ s \r,d)\m) = E P^ n ~ m) f s nm {r) ■ 

m,n=0 m,n=Q 



2 The function (s,z) i— > lj(s,z) :— tr[pG' s ^ (z)] is also known as the photon number tomogram [20l I21j. that 

Is) 

is, the tomogram is identified with the collection of probability densities G p , s G N. 



Using equation (DO), the explicit form of the functions /* m can be written as 



f s / x r _ n max{0, S -n}+max{0, S -m} / mhl{n, S }\ mhl{m, S }\ 

/ "" ' W max{n,s}!max{m,s}! 



-r 2 | S -n| + |s-m| r [s-n| / 2\ r [«-m| / 



The mapping 6* i— > G^ s '(r,9) is weakly continuous for each r G [0, oo), and ||Gl^(r, 9)\\ = 1 
for all r G [0, oo), 9 G [0, 2tt), so the operator 

Gf(r) := — 1^ e iW G ls) {r,9)d9 



2vr ... 

is well-defined as a weak integral. In addition, we have 

G[2(r) :=tr[pGi s) (r)] = -L e* w G^(r,9) d9, 
for all states p. A simple calculation gives us 

oo 

G l ;l(r) = Y,Pn + iAMD(r)\s)(s\D(r)*\n + l), 

n=0 

for all r G [0, oo). 

The probability distributions G$ , also known as displaced photon distributions, are closely 
related to the A-parametrized phase space quasiprobability distributions, first presented by 
Cahill and Glauber [TIE]- To clarify this, let us recall the definition of these distributions. For 
each A G C, |A| < 1, define the operator density W x : [0, oo) x [0, 2ir) — > C(H) by 

oo 

W\r, 9) := (1 - A) J2 ^ k D(re ie )\k)(k\D(re l9 y, 

k=0 

and the corresponding probability density W x by Wp(r,9) = tr[pW x (r, 9)). It is clear from 
these definitions, that indeed, one has 

oo 

W p \r,9) = (1 - A)^A fc Gf(r,0). 

To obtain the displaced photon distributions from the A-distribution, we first note that 

|G'*>(r,0)| = \{k\D(re»y P D(re»)\k)\ < \\p\\ = ||p||i < 1, 

for all r G [0,oo), 9 G [0,27r), which follows from the Cauchy-Schwarz inequality. This then 
implies that (1 — A)~ 1 W / ^(r, 9) = Ylk^o A fc Gp^(r, 9) is a power series with respect to A, converging 
absolutely for all A G C, |A| < 1, suggesting that the series can be differentiated around the 
origin term by term. A direct calculation now gives us 



A=0 ' u—n A=0 



Q.s X k 



A=0 fc=0 



s\d\ s 
since -gxr| A=0 = sld sk . 

These of course give us, at least in principle, the possibility of constructing either of the dis- 
tributions from the other. In a recent paper [14], rigorous proofs for two reconstruction formulas 
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for the A-distributions were given. In view of this, the knowledge of all of the distributions G p , , 
s G N, allows state reconstruction via a detour. 

4.1. Measuring the displaced photon distributions. We will now review the possibility 
of measuring the G'^-distributions with an eight-port homodyne detection scheme. For a basic 
reference concerning the setup, see e.g. [H]. In [13] a rigorous proof was given for the fact that 
with this scheme, any covariant phase space observable can be obtained as a high amplitude 
limit. The detector consists of two pairs of photon detectors and the amplitude-scaled photon 
differences D\ and D2 are measured. Four input modes are involved; the signal mode, a vacuum 
mode, a local oscillator in a coherent state, and a parameter mode which defines the observable 
to be measured. If the parameter mode is in a state S, then the phase space observable E csc , 
where C is the conjugation map ip > (x t— > if)(x)), can be obtained as the high amplitude limit 
(see [13]). 

The first obvious way to measure the G'^-distributions is the direct measurement in the 
sense of the above limit. However, this requires ideal detectors, and the parameter field needs 
to be prepared in a number state The preparation of the number state is highly untrivial 

and is by itself an active area of research. Several theoretical models, mostly in the context of 
cavity quantum electrodynamics, for the preparation of an arbitrary number state have been 
proposed (see e.g. [HI EH [3], [5])- Even though this gives a theoretical method for measuring 
the distributions, it is not a practical one since the preparation procedures work only for small 
photon numbers. To avoid the problem of number state preparation, we consider an alternative 
point of view. 

Consider the measurement of the Q-function of the electromagnetic field by means of the 
above experimental setup. In this case the parameter field is in the vacuum state |0)(0|. If the 
detectors are non-ideal, with a detection efficiency 77 each, the measured distribution is actually 
the A-parameterized distribution, with A = 1 — 77 [T7l [9]. Suppose now that the detector 
efficiencies are close to unity, that is A ~ 0. Then, by adding suitable beam-splitters into the 
measurement scheme, one is able to measure the distributions corresponding to the parameter 
A', for which A' > A. An equivalent scheme would be one where the detector efficiencies could 
be adjusted. Proceeding in this manner, one obtains a function A i— > W^. In an ideal situation 
where 77 = 1 one could thus differentiate this s times with respect to A, and obtain the G^- 
distribution according to equation ([!]). Even in the nonideal case, one can obtain some kind of 
an approximation for the G'^-distributions, provided that the A-dependence of is regular 
enough to allow an extrapolation to the values close to the origin. 

4.2. Reconstruction from the set {Gp \s 6 N} of distributions. If one has knowledge of 

\s) 

all of the distribution G p , s E M, recovering the diagonal elements of the density matrix is a 
trivial task. Indeed, by definition one has 

GW(0) = tr[pD(0)\s)(s\D(0y] = (s\p\s) = p ss , 

suggesting that in order to reconstruct the diagonal elements of the state matrix, one needs 
to measure the observable around the origin for all s £ N. The reconstruction of the 
off-diagonal elements is a more complicated matter. 
Let I 6 N, I > 1, so that 

00 

° l pl( r ) = ^2pn+l,nfn, n+ l(r), 
n=0 

10 



where the functions fn n+ i were defined in equation ([3j). Define a function g% : (0, oo) — > C via 
9i{f) = e r V~'. Suppose that I < s. Then the limit T l sn := \im r ^ Q (gif° n+l )(r) exists, because 
\s — n\ + \s — (n + l)\ > I for any n G N, s > I, and can easily be computed using the fact that 
^(0) = ( m + Q ); the result is 

0, n < s — I; 

T ln = \ (-^y/W is-nyln+l-sr ^ ~ I < U < S; 

, n > s, 

In addition, assuming n > s and r G (0, 1), and using the fact that |L™~ s (r 2 )| < (™)e2 r2 [TJ p. 
786, 22.14.12] we get 

I / \ rs 1 M / e / nS (n + iy 

\9l(r)fn,n+l( r )\ < 



s! V (n — s)\ (n + I — s) 



which goes to zero, as n — > oo. This implies that sup raeNrg( - \gi(r)f^ n+l (r)\ < oo. Since 
Yln=o \Pn+i,n\ < 1, it follows that the series £)^1 p n +i, n gi(r)f n ,n+i(r) converges absolutely and 
uniformly on the interval (0, 1). Thus, the limit 

oo 

]imgi(r)G^(r) = lim p n +l,n9l(r)fn , n +A r ) 

n=0 

may be taken termwise. This gives us the infinite matrix identity 

oo 

d l s ■= lirn \gi(r)G l ^(r) = Vrj n p n+ijn , 

n=0 

which, in this case, holds for all states p. Inserting the explicit form of T l sn we obtain 



4 = E (-D-V^ 



1 



n! (s — n)\(n + Z — s)! 

n=s—l 



s+l / ^7] j 



Defining c* n := \f^mPn,n-i for n > / and c* n = otherwise, and a l sn := (-l) n *(JJ. 



(n-QI 

we can write equation (jHj) as 



^ ] a sn C n-> 
n=0 

since s > I. The infinite matrix (a^ n ) 8ifie N is now of the type considered in lemma [3] with 
a l u = (— 1) M (^). Consider now the sequence (b l u ) ue ^ with b l u = ("t^ 1 )? and the infinite matrix 
(^ns)n,seN) = ^s-n- This is an upper triangular matrix, and we have 

°° °° /A // - 1\ 

^ ^ Q'srfins ^ ^ ^ns^sn ^nrfinn ( 0/ V i 1/ ^' 

n— n o— n \ / \ / 



n=0 s=0 
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To calculate the off-diagonal elements of the product matrices, let k > 1. Using formula (5) on 
page 8 of [23], we find that 

= E(-^(„i.)C +i 7- B i + ' _1 ) 

n=0 n=s \ / \ / 

-d-^OP-^- 1 )-?; 1 )- 

For the other case we use formula (5d) on page 10 of [23] to obtain 

Since the lower diagonal elements of the product matrices are zero by definition, we find that 
( a sn)s,neN an d (6[ ls ) n ,seN are formal inverses of each other. 
Suppose now that 



(6) lim ma p m , m -i = 

m— >cx) 

for all I G N. Then the limit condition of lemma [3] is satisfied for each n G N, since 



/ / j- / m\ (m + l — k — 1)1. 

\bm~kCm\ = /TV (rn-l)\ {l-l)\{m-k)V P ' 



m,m—l 



< j^j - ^ m?(m + I - k - l) 1 l \p m ,m-i\ 

< 2^ 1 mi / ~ 1 |p mjm _ i | 

for m > I. Under this condition, we then have the convergence 

oo oo 
s=0 s=n 

or equivalently 

c n+l — 

s=n+l 

which gives us the reconstruction formula 



<7) ^=(-d'«^e +i ( s 7-"7>-' 

where d! s is a quantity which can be calculated directly from the measurement statistics. 

Remark 2. Notice that the condition (JHJ) for a given I G N is a sufficient condition for the recon- 
struction of the Ith diagonal of the density matrix. For I = 0, for example, the reconstruction 
formula works for all states p. However, in the general case, the validity of the formula depends 
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on the state in question. To illustrate this fact, consider the vector state ip = -j= Y^=i „l n )- 
Condition (El) now states that we should have 



lim 



777,2' 



0. 



n-»oo 777(777 — I) 

which is clearly not true for I > 2. It is easy to check that even the weaker condition, namely the 
limit condition of lemma [3l is unsatisfied. This then suggests that the reconstruction formula 
(|Zj) does not work for all states. 



I«> 
p 1 



4.3. Reconstruction from a single distribution. If we want to use a single distribution G 
the reconstruction formula becomes more complicated, and we were not able to satisfactorily 
solve the convergence issues in the case of an infinite density matrix. Consequently, we will 
assume in the sequel that the matrix is finite. This corresponds to the discussion of section [3] 
concerning finite sequences. 

The reconstruction makes use of the connection to the A-parameterized distributions W*. It 



follows from equation 



, that 



1 d s 



((i-A)-X*M)) 



A=0 



for all I G N. On the other hand, we have for all states p and I G N 



1 r 2w 00 

W£,(r) := - / e* w tr[pW\r, 9)] d9 = J2 Pn + i, n K* n 



+i\T), 



where 
(8) 



K, n+l {r) 



1 - A)^A fc (n|D(r)|A;)(A;| J D(r)> + /). 



k=0 



This series can again be differentiated termwise, and we get 



n=0 



1 d s 



((1 - A)-^ n+ ,(r)) 



A=0 



The explicit form of the functions K ntTl+ i is given by the formula of Cahill and Glauber [7] 

-(1 - A) i+1 e- (1 - A)r V'A n L^((2 - A - A^r 2 ) 



K,n+i(r) 



77! 



(77 + /)! 

n 



u=0 



-1 \\2u+l+l\n—u„2u+l 

1 ~ A ) A r C ~(l-X)r 2 

(n - u)\(l + u)\u\ 



Before proceeding any further, we prove the following lemma. 

Lemma 4. Let k, p, q, s E N and xGl 

(a) 

1 d k {l - A) p A 9 



fc! 



d\ k 



1 



k+q 



A=0 



P 

k — q 



which is if and only if k < q or k > p + q. 
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(b) 



1 d s (l - \) p \ q e 



Xx 



Si 



dX s 



wm{s,p+q} , _^ k+q 



A=0 



E 

k=q 



(s — k)\ \k — q 



P 



which is for all x if and only if s < q. 
Proof. By direct calculation we get 



d k (l - X) p X q 



d\ k 



E 



A=0 t =Q 



k\ d l X q 



t dX l 



d k ~\i - xy 



A=0. 
qld q ,t 



dX k ~ 



A=0 



from which (a) follows. Part (b) follows from (a) and the calculation 



1 d s {\ - X) p X q e 



Xx 



SI 



dX s 



A=0 



. /s\ d k (i - xyx q 



k=0 



dX k 



A=0 



d s-k e Xx 

dx s - k 



A=0 



2^ ( s - 



k=0 



Uk+q 

s — k)\\k — q 



Now fix s E N and denote x = r 2 . We have two different cases depending on whether I 
even or odd. We will start with the even case. 

The even case. Suppose that I = 2h for some h EN. Then, by lemma 0] we get 



1 d s 
^.dX? 



((1 " AJ-^^Ca/E)) ^ n = e-^n\{n + 2h)\ £ 

A — 



it=max{0,n- s} 



u\(n - u)\(u + 2h)\ 



x 



x 



mm{s g 2/l+n} (-l) fc +-V 2(u + h) 

k=n—u 



(s — k)\ \k — (n — u) 



X 



u+h+s—k 



For any t E N, define 



A=0 



i=0 



so that 



x=0 



n=0 
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Now 



^ y/n\(n + 2h)\ min ^ 2fc+ "> ^ 2(« + ft) \ frx'** 1 *-* 



u!(n — u)!fu + 2/i)! ' (s — fc)! \k — (n — u) J dx l 

= t! <5t,u+h + s-fe 



x=0 



^ t! v /n!(n + 2/i)! (_i)^+*+n ^ 2(u + ft) 



u!(n - u)\(u + 2h)\ (t-h-u)\ \h + s-t-n + 2u 

«=max{0,n-s} v ' y ' v ' v 

t! v /n!(n + 2/i)!(-l) s+h+ * + ™ min g- /l} i (2( !i +/i)Vs-n + u 



f/i + t + n — s)\ ' (n — «)!(« — n + s)! \ tt J \t — h — u 

: ; /u=max{0,n-s} v 7 v y x 7 x 



: iff h+t+n-s<0 v 

= Oifft-k<Oori-k<n-sorn<f-/i-s 



Since H 2h (t, n) = if t < h we next assume that t > h and get 



at s-M-t 



(9) ^G^{y/x) x=q = H s 2h {t,n)p n+2h , n . 

n=max{0,s- h— t,t— h— s} 



Let t > s + h and denote p = t — s — /i. We have 



e^Sjv^) = J2 H s 2h (p + s + h, n)p n+2h>n . 

n=p 



Define an upper triangular matrix {A s '^ h ) p tn£ ^ by 



A s p n h := H° h (p + s + h,n), n > p. 
According to Corollary it has an inverse matrix 

/ TDS,2h\ 

\ U np )n,p£N- 

Using this, we get the state reconstruction formula 

00 PjP+s+h 

(10) Pn + 2 h ,n = £ BS 4 h ^T h ^U^) 

p=0 

The odd case. If I is odd, that is, I — 2h + 1 we get 

1 d s 



x=0 



s\d\ s 



((1 - A) ~ 1 K^ n+2h+1 ( yfx ) ) = V£e-Vn!(n + 2/i + l)! 



A=0 



1 ^ {s ^ +n} { _ 1)k+n _ u f2{u + h) + ^ 

2-^1 ii\(r, — iA\(n -I- Oh 4- "HI 



_ { o,n- s} < n ~ U ^ U + 2h + ^ (s-W \k~(n- u) 

For all t G N, define 



x 
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A=0 



x=0 



and calculate 



t\^n\(n + 2h + 
^ +l(t ' n) " {h + t + n-s + l)\ X 



= iff ft+t+n-8+l<0 
min{n,t— h} 



u=max{0,n- s} 



(n — u)\(u — n + s)\ 



2(u + h) + 1\ /s - n + u 



u 



t — h — u 



= Oiffi-/i<Oorf™/i<n-sorTi<(-/!-s 

Since H% h+1 (t, n) = if t < h we next assume that t > h and get the infinite matrix identity 

^e*G% h+x {Sc) 



(11) 



x=0 



22 H2 h+1 (t,n)p n+2 h+l,n- 
n=voax{0,s—h—t—l,t—h—s} 

Similarly to the even case, assume t > s + h and denote p = t — s — h to get 



Qp+s+h 



dx- 



2s+p 



x=0 



Yl H 2h+i(P + s + h,n)p 

n+2h+l,n- 



n=p 



Let (B^p +1 ) n , P eN be the inverse of an upper triangular matrix (Ap^ h+1 ) Ptn( zfq with 



A% h+1 :=H°(p + s + h,n), n>p. 



Thus, we get the formula 
(12) 



°° 2h+l d p+s+h . \ a \ 

Pn+2h+l,n = Y Kp h+1 Q x p+s+h G p,2h+l ( V^) 



p=0 



x=0 



Example 1. As an illustrative example, we consider the simple case of the observable 
generated by the first number state |1)(1|, and the system in the state p = Ym:=o a k\k){k\, 
where a*. > for k — 0, 1, 2 and Y^1=o a k — 1- The density matrix is now diagonal, so we use 
the reconstruction formula (fTUl) for s — 1, h — and n — 0, 1, 2. We easily obtain the quantity 



e x G l pl(x) = a x + a x {l - x) 2 + \a 2 x{2 - x) 2 , 



so that the nonzero derivatives are 



«o — 2«i + 2a>2, 



x=0 



—e x G ll) (x) 



2a± — Aa 2 , 



d 3 



e x G l ! ] Jx 



x=0 



dx 3 pfi 



3a?. 



x=0 



The reconstruction formula can then be written explicitly as 

(13) pnn = (ao - 2«i + 2a 2 )B^ + (2a>i - Aa 2 )B l n ® + 3a 2 ^2°- 

To obtain the matrix elements we first calculate the matrix (Ap^) p neN . It follows from 
elementary calculations that 



A 



1,0 

pp 



p+1, A^ p+1 = -(2p + 2), 
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A 



1,0 

P,P+2 



P + 2, 



and 



otherwise. In matrix form, this reads 



/ 1 


-2 


2 














2 


-4 


3 








n 

W 


n 






4 


n 











4 


-8 


5 














5 


-10 

















6 















The inverse matrix can now be calculated, and we obtain 



[B 



1,Q\ 
np ) 



( 1 1 | J 

I 2 I 

u 2 3 8 
| i 





Villi 



J5 



2 ^ 5 



_31 

w 

1 
2 1 4 



5 





/ 



Inserting the appropriate elements into equation (1131 ). we find that 

, 2 



Poo 


= 1 • 


(«0 


- 2ai 


+ 2a 2 ) 


+ 1- 


(2ai 




= 0- 


(«0 


- 2ai 


+ 2a 2 ) 


1 

+ 2 


•(2a! 


P22 


= 0- 


(«0 


- 2ai 


+ 2a 2 ) 


+ 0- 


(2ai 



■ 3a 2 = «o, 



4a 9 



i-3a 2 : 



«2, 



so that the formula does indeed give the correct values for the matrix elements. 

5. Discussion 

In this article we have considered the problem of reconstruction the unknown state of a 
quantum system from the measurement statistics of phase space observables generated by the 
number states. The two reconstruction formulas we have derived are dealing with two different 
scenarios and, as such, have obviously different advantages and disadvantages. 

In the first case, the reconstruction requires measurements of all of the observables . From 
the practical point of view, this is of course impossible for many reasons. First of all, in the 
eight-port homodyne detection scheme, for the measurement of E' s ', a parameter field needs to 
be prepared in the number state |s)(s|. At the present, this is possible only for small values of 
s. Nevertheless, it might be reasonable to expect that future progress could allow sufficiently 
large number state preparations, so that the reconstruction would be possible with adequate 
precision. An advantage of this method is that one only needs to measure the observables near 
the origin, and the whole phase space does not need to be scanned. 

The second method uses a single obervable and thus the problem of number state gener- 
ation for arbitraily large values of s is removed. This time the region of phase space needed is 
considerably larger, since, in principle, the method involves integrals over the whole space. This 
is of course an idealization which can not be considered in practice. From the numerical point 
of view, this method seems to be very manageable. Suppose that one measures the observable 
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for some s. Using polar coordinates for the resulting phase space ditribution, one can then 
integrate the distribution with respect to the angle variable over [0, 2n) to obtain the quantity 
gJ^q. If the density matrix is assumed to be finite, say, an N x iV-matrix, then e r2 G^ Q (r) is a 
polynomial of order 2N + 2s — 2. Choosing a sufficiently large N E N, one can thus fit such 
a polynomial into e^G %{r). This also fixes the size of the approximative density matrix. In 
view of the reconstruction formula, this assumption of finiteness has the crucial consequence 
that only the inverses of iV finite matrices are needed. Similarly, the quantities G^\ are needed 
only for 1 = 0,...,N- 1. 
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